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We study a quasi-two-dimensional electrostatic drift kinetic system as a model for near-marginal ion tem¬ 
perature gradient (ITG) driven turbulence. A proof is given of the nonlinear stability of this system under 
conditions of linear stability. This proof is achieved using a transformation that diagonalizes the linear dy¬ 
namics and also commutes with nonlinear E x B advection. For the case when linear instability is present, 
a corollary is found that forbids nonlinear energy transfer between appropriately defined sets of stable and 
unstable modes. It is speculated that this may explain the preservation of linear eigenmodes in nonlinear 
gyrokinetic simulations. Based on this property, a dimensionally reduced (oo x oo —>• 1) system is derived 
that may be useful for understanding dynamics around the critical gradient of Dimits. 


a. Motivation. The most unstable ion temperature gradient (ITG) modes in toroidal fusion devices, such as 
tokamaks and stellarators, rely on a drift resonance due to curvature in the magnetic guide field. One signature of this is 
a concentration of turbulent activity around regions of “bad curvature.” This reflects the linear instabilities themselves 
which peak or “balloon” in these regions. For sufficiently unstable modes, two conditions can be simultaneously melii: 
(1) the ballooning is strong (the mode is strongly localized) and (2) the timescale associated with ion propagation 
along this structure is long compared to the mode frequency {i.e. the parallel transit frequency is small), rendering 
the dynamics nearly two-dimensional. For these modes, the scale lengths perpendicular to the magnetic field are 
significantly larger than the ion Larmor radius so finite Larmor radius (FLR) effects can be considered small. 

If all of the conditions above are met, even if only for a range of spatial scales, it is appropriate to study the two- 
dimensional drift kinetic equation to understand the dynamics that results. Furthermore, a precise understanding of 
the basic properties of this equation is, generally speaking, essential as a basis for the understanding of the general 
gyrokinetic system, of which it is a limiting case. 

There are two important fundamental questions that further motivate this work. First, why does plasma turbu¬ 
lence sometimes retain the features of the linear instabilities that drive it? In the fully nonlinear state, one might 
expect the plasma fluctuations to become disordered to the extent that linear modes cannot be recognized by their 
spatial structure, frequency and growth rate. Yet, all of these features have indeed been observed in fully developed 
turbulence^^— , even when it is strongly driven^. Second, what determines whether turbulence will exist in the absence 
of linear instabilities? Whenever there is non-uniformity of the background, there is a source of free energy to excite 
turbulence, and many plasma systems are known to support turbulence even when there are no linearly unstable 
modes to drive iti^^— It should therefore be valuable to identify linearly stable cases in which turbulence provably does 
not exist, and to understand the underlying reasons. 

b. Basic analysis. Let us consider the following two-dimensional drift-kinetic equation for the ion species of a 
magnetized plasma (obtained from the gyrokinetic equation by taking the Larmor radius to be small and ignoring 
parallel streaming): 


dt 


+ ikyVdf + iky{vd 


- 'y*)¥>Fo = ^ekpqV3*(p)/*(q), 

p.q 


(1) 


where Lp = qifji/Tc \s the normalized electrostatic potential, Fq is the background ion distribution, and the perturbed 
distribution function is f{'k.,v±,v\\,t) where k = xfca, -I- yky. The nonlinear coupling coefficient is eicpq = <5(k + p + 
q)To/((3'ii?)z • (p X q) is antisymmetric under the exchange of any two components, i.e. has the symmetry of the 
Levi-Civita tensor. The electrostatic potential is determined by the quasi-neutrality constraint, assuming a modified 
Boltzmann electron response. The quasi-neutrality constraint is written 


q noo 1-00 

( 1 - 4 j ,) t (/5 = — / dv\\ Vj_dv_i_f, ( 2 ) 

^0 J — oo J 0 

where 5k is the discrete delta function. Note that this implies the zonal ion density is zero (at this order in the small 
Larmor radius expansion). To close this system we need an equation for the zonal electrostatic potential [ky = 0), 
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which is given by the vorticity equation (obtained from the gyrokinetic system, including terms of order k^p^; see 
Plunk, Tatsuno, and Borland— Eqn. D.9) 

ekpq<p*(p) ['?V*(q) - p • q 7 "I(q)], ky = o. ( 3 ) 

p,q 


We use the following additional conventions: Fq = no(wT^7r) exp[— v'^ = + z;jj, r = Tio/{ZTeo), 

Z = qi/\qe\, To = Tio, T± = (Tono)“^27r/ 'y_Ld'y_L(iu||/(mt>]^/2), Vd = p\/2/(yTR)[v\l2+v‘^^ and u* = pv^/{LnV2){l + 
p[v^/ vt^ — 3/2]), L~^ = din no/da;, = dlnTg/dx, 77 = Ln/Lx ■ 

In the presence of a sufficiently large temperature gradient and magnetic curvature drift of the correct sign (i.e. 
bad curvature), this system is linearly unstable to ITG modes of the “toroidal branch.” A stability analysis of this 
mode was done by Biglari, Diamond, and Rosenbluth— who found a sufficient condition for stability, Lx/Ln >3/2. 
Although a necessary and sufficient stability criterion could not be obtained analytically, the quantity k = R/Lx is 
typically cited as the important instability parameter governing the toroidal ITG branch. For the present purposes, 
we need not specify the value of the parameters in our system; we will simply specify whether or not the linear ITG 
mode is stable. 

Now there is a simple quadratic form that is a dynamical invariant of the system [T][21 It is found by multiplying 
Eqn. [ 1 ] by the complex conjugate /*(k) and then taking the real part of the result; this annihilates the Vdf term. 
Then, we divide by the function «:(7 ;_l,?;||) = {vd — 'i’*)Fb, integrate over velocity, and use Eqn. [5] to annihilate the 
term that goes as ikyf*Lp + c.c.. Finally, we sum over k, and use the antisymmetry property of the nonlinearity to 
show that its total contribution is zero. The result is 


dt 


27r / v±dv±dvi\ 

1 . 


\f? 


= 0 . 


(4) 


This is a useful constraint if the function n is sign definite, because then the integral defines a norm for the function 
/. In this case Eqn. Sj is a statement of nonlinear stability. The sign-definiteness of k can be guaranteed if the 
stability criterion LxjLn > 3/2 is satisfied, but this is not a very useful criterion as already mentioned. In what 
follows nonlinear stability is proved in a much stronger sense, i.e. by only assuming linearly stability it is shown that 
Eqn. [ 1 ] is also nonlinearly stable in the sense that an appropriately defined norm is invariant in time. 

c. Nonlinear stability. We now prove a nonlinear stability theorem by explicit use of the transformation that 
diagonalizes the linear operator. The proof relies crucially on the existence of a complete set of eigenmodes of the 
linear operator, which we now assume to be stable. Because the dependence of this operator on k is only in the 
overall multiplicative factor of fcy, it should be immediately clear that the eigenmode analysis is independent of k. 
Gonsequently, linear diagonalization commutes with the nonlinear E x B advection and the transformed system is 
manifestly stable. Let us show this in detail. 

First, we reduce the velocity space of system [T] by using the change of variables (n_L,?;||) —>• (il,A). The idea is 
that the linear phase-mixing in Eqn. [T] is really one dimensional, so we should be able to eliminate the ignorable 
dimension. The transformation can be done in two steps, defining first the polar coordinates w = {w'j_ + and 

9 and corresponding cartesian coordinates w± = v±ly/2 and W|| = z;||, such that 6 is the angle between the r(;_L-axis 
and the vector (wx, wn). Then we define v = Vd = jcd and A = sin(0), where Cd = vxR/{^/2p). The equation may 
be integrated over A to obtain the reduced system in ^(k, v, t) = dA/: 


dg_ 

dt 


-I- ikyvg + ikyG{v)n = (p)g*(q), 

Pq 


(5) 


where 


n = 


dvg{v), 


( 6 ) 


and 
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k{v, a) 


exp 


Cdv{2 - A^)'' 


tiT 
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In the new variables we write k(z;, A) = v — c^,{l + r]{{2 — X^)vcd/vT^ — 3/2))] with c* = pvTI{LnV2). The linearized 
system has a continuum of neutrally stable solutions, and in the present form it is almost directly amenable to 
solution by the Morrison transform j^^d^ see also Plunk— for a recent application of this technique. However Eqn. [5| 
is different in form from the classical problem of Landau damping in that the velocity-space domain is only semi¬ 
infinite, V G [0,oo). The continuous spectrum is also only semi-infinite, and the required mathematical tools (Hilbert 
and Fourier transforms) are defined on an infinite domain. To circumvent this difficulty, we can simply extend 
the space v to the full domain. Formally we consider the new system on u G (— 00 , 00 ) for the function 
n^^\v) = J^^dvg^^\v) and = G(u)0(u), where 0 is a step function, 0[u] = 1 for u > 0 and otherwise zero. 

The linear dispersion relation for the discrete spectrum is equivalent to that ofO and the nonlinear solutions of the 
original system are a subset of the solutions of the extended system. In what follows we will omit the superscripts of 
the extended system to keep the notation simple. 

Let us define the Morrison transform in the notation of Plunk— 



9- 


( 8 ) 


where D± = 1 ± 2TriG±. The positive- and negative-frequency parts of an arbitrary function g{v) are defined in terms 
of the Fourier transform by 


^±00 

g±{u)=± dve^''^ 

^0 


dv- 


27r 




The Morrison transform can be directly applied to Fqn [5| to obtain 


(9) 


-f ikyug = ^ekpqV?*(p)5*(q). (10) 

p,q 

The following conservation law is then obtained, completing our proof of nonlinear stability: 

k 

Note that this applies point-by-point in velocity space it, in contrast to the quantity in Fqn. |4| Fqn. [TT] constitutes 
a statement of stability because Fqn. |8| is an invertible transformation. Invertibility is equivalent to completeness of 
the eigenmodes, which, as explained by Kampen and Felderhof— , is guaranteed if the function £)+(it) has no zeros 
in the upper half plane. This is true if and only if the dispersion relation for the linear system has no roots in the 
upper half plane. 

d. Unstable modes and “energetic isolation.” It turns out that the above analysis may be extended to 
cases with a discrete spectrum, resulting in a theorem about energy flow in the unstable system. Following the work 
of Case—, we may use the orthogonality of the eigenmodes to separate the continuous and discrete parts of the 
eigenspectrum. Numbering the discrete modes 1,2,3,..., m, we define the projection operator for each mode 


'Pi [ 9 ] = j dvh,{v)g{v), ( 12 ) 

where gi{u) = —G{u)/{u — Ui) is an eigenmode of Fqn. [5l hi{v) = —l/(u — Ui) is the corresponding eigenmode of an 
adjoint equation, and Ci = dvG{v)/{v — Ui)^; see Case— for details. Note that the discrete spectrum generally 
includes pairs of damped and growing modes, and also stable non-singular eigenmodes (exceptional cases). The 
projection operator may be applied to Fqn. [S] to obtain a system that describes the linear evolution and nonlinear 
interaction among the sub-population of modes {^^(ki), 5 i(k 2 ),...}, indexed by the wavenumbers of the system. That 
is, writing ^(k, v) = f A{u, \i)g^(v)dv -\- O‘i0^)9i{^) have 
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9aj(k) 

dt 


+ ikyUiai{k) = y^ekpqy’*(pX(q), 
p,q 


( 13 ) 


implying 


^ (14) 

k k 

We may also isolate the continuous spectrum by simply applying the operator V' = 1 — 1° Eqn. [5j This 

yields an equation for g' = 'P'{g\ (of the same form of Eqn. [5]) that only has a continuous stable linear spectrum. To 
this system the invertible Morrison transform may be applied as done in the previous section. The sum of the energy 
associated with these modes is an exact invariant of the full nonlinear equation, i.e. Eqn. [TT] for g'. In summary, 
we may decompose the system into sets of modes (i.e. the sets {gi^ki), gi{'k 2 ) ,...} and {g,^(ki), g^(k 2 ), ...}) that can 
be dynamically coupled via the nonlinearity but are energetically de-coupled. Note that for a set of unstable modes, 
the zonal components ai{kx,ky = 0) are linearly stable, but can be excited nonlinearly. This energetic isolation of 
linear mode populations is an intuitive property: energy is passed between equivalent modes because (1) the linear 
eigenmodes do not depend on the wavenumber k and (2) nonlinear mode coupling preserves the velocity dependence 
of g. 

Note that a consequence of energetic isolation is that damped eigenmodes^ cannot be excited by nonlinear interaction 
with the set of unstable eigenmodes. This does not contradict studies of fluid and gyrokinetic system a^^d'^ that 
demonstrate a significant degree of such excitation, because those systems have additional mechanisms, e.g. nonlinear 
phase mixing, that enable energetic exchange between different mode types. However, the present work does imply 
that such energy flows should be small to a degree determined by how well the assumptions of our system (namely 
k'^p^ <C 1 and small collisional and parallel transit frequencies) are satisfied. 

e. Dimensional reduction The exact energetic isolation of eigenmodes suggests that a dramatic simplihcation 
could be attempted. Let us assume the stable modes can be neglected (in a real plasma, an additional mechanisms 
like collisional dissipation could drain the energy of the stable spectrum) and that there is precisely one unstable mode 
for each wavenumber k with ky ^ 0, which, dropping the indices, we notate as a. Then by Eqn. [2l the non-zonal 
part of the electrostatic potential is proportional to the amplitude of the unstable mode, i.e. ^/^(k) = 27r(cd/T)a(k) for 
ky 7 ^ 0. Now we can divide the nonlinear coupling terms on the right hand side of Eqn. [13] into three groups, (1) the 
non-zonal/non-zonal interaction terms, Py 0 and qy ^ 0, (2) interactions involving the zonal potential, py = 0 and 
Qy 7 ^ 0, and (3) interactions involving the zonal part of a, Py ^ 0 and qy = 0. The first group of interactions clearly 
sum to zero by the antisymmetry of e. The final group must also be zero because a(k) = 0 for ky = 0, which follows 
from Eqn. [2] and the fact that the eigenmodes have unit density. What remains is a very simple system consisting of 
two equations. The first is 


5a (k) 
dt 


ikyUiaik) 


X! ekpq:/J*(p)a*(q), 
T^y —® ^Qy 


(15) 


where the sum is performed over the set of p and q satisfying py = 0 and qy ^ Q . This equation is then closed by 
Eqn. 13] for the zonal vorticity, which can now be written in terms of the drift wave amplitudes: 


= X!^kpqa*(p)a*(q) - a/3p • q] , fcj, = 0. (16) 

p,q 

where a = 21:04 It and 

P = [ x^j_dx± [ dx\\ (— —exp(—— Xm). (17) 

Jo J-00 \Vd-u^J " 

Eqn. [15] describes the nonlinear zonal shearing of a set of unstable drift waves, and the voriticity equation [16] 
describes the nonlinear response of the zonal flows to the bath of drift-waves. Note that some kind of dissipation is 
needed to make this system well-behaved since there is no mechanism to saturate the unstable modes: if there are 
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finite amplitude drift waves a(k), the set will always grow in energy. If small-scale dissipation were to be added, 
then zonal shearing would be able to transfer energy to large sink and presumably cause saturation. Note that 
the non-density components of the zonal modes become larger in amplitude, and dynamically important above the 
nonlinear critical gradient, and are responsible for nonlinear zonal flow decay via the tertiary instability^^. Thus it 
seems likely that the condition a(k) = 0 for ky = 0, restricts the regime of validity of this system to around or below 
the nonlinear critical gradient of Dimitai^. 

/. Conclusion. We have shown that the two-dimensional drift kinetic system, Eqns.[T][31 satisfies two properties: 

(1) it is nonlinearly stable when linearly stable and (2) it generally only allows equivalent linear eigenmodes to exchange 
energy nonlinearly. The first property means that subcritical turbulence cannot be excited unless by a mechanism 
not present in our system, e.g. by three-dimensional mode coupling^ or the parallel velocity gradient drived. This is 
an encouraging fact because it seems to limit the contribution of the magnetic drift resonance as a turbulent drive to 
only cases where it can induce a linearly unstable mode. The second property, i.e. the fact that eigenmode sets are 
energetically isolated, may help explain observations^^ that strongly-driven ITG turbulence preferentially dissipates 
energy by cascade to smaller scales, rather than by nonlinear transfer from an unstable mode to a damped mode 
at a similar scale. Indeed, as the turbulence is more strongly driven, it peaks at larger scales, kp 1, making our 
drift kinetic system more applicable. Energetic isolation may also help explain why gyrokinetic turbulence retains 
recognizable features of the linear eigenmodes, even when strongly driven. The preservation of linear modes is a 
mysterious but well-known property^r— of gyrokinetic turbulence with no accepted explanation. Crucially, we have 
not assume weakness of the nonlinearity, so this result constitutes a new kind of mechanism for preserving linearity 
that might motivate quasilinear models when the assumption of weak nonlinearity is unjustified. 
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